home *** CD-ROM | disk | FTP | other *** search
/ STraTOS 1997 April & May / STraTOS 1 - 1997 April & May.iso / CD01 / INTERNET / SITES / LITTLE / P3SRC.ZIP / ATARI / RADIOSIT.C < prev    next >
Encoding:
C/C++ Source or Header  |  1997-01-17  |  31.8 KB  |  1,048 lines
  1. /****************************************************************************
  2. *                   radiosit.c
  3. *
  4. *  This module contains all radiosity calculation functions.
  5. *
  6. *  This file was written by Jim McElhiney.
  7. *
  8. *  from Persistence of Vision(tm) Ray Tracer
  9. *  Copyright 1996 Persistence of Vision Team
  10. *---------------------------------------------------------------------------
  11. *  NOTICE: This source code file is provided so that users may experiment
  12. *  with enhancements to POV-Ray and to port the software to platforms other
  13. *  than those supported by the POV-Ray Team.  There are strict rules under
  14. *  which you are permitted to use this file.  The rules are in the file
  15. *  named POVLEGAL.DOC which should be distributed with this file. If
  16. *  POVLEGAL.DOC is not available or for more info please contact the POV-Ray
  17. *  Team Coordinator by leaving a message in CompuServe's Graphics Developer's
  18. *  Forum.  The latest version of POV-Ray may be found there as well.
  19. *
  20. * This program is based on the popular DKB raytracer version 2.12.
  21. * DKBTrace was originally written by David K. Buck.
  22. * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
  23. *
  24. *****************************************************************************/
  25.  
  26. /************************************************************************
  27. *  Radiosity calculation routies.
  28. *
  29. *  (This does not work the way that most radiosity programs do, but it accomplishes
  30. *  the diffuse interreflection integral the hard way and produces similar results. It
  31. *  is called radiosity here to avoid confusion with ambient and diffuse, which
  32. *  already have well established meanings within POV).
  33. *  Inspired by the paper "A Ray Tracing Solution for Diffuse Interreflection"
  34. *  by Ward, Rubinstein, and Clear, in Siggraph '88 proceedings.
  35. *
  36. *  Basic Idea:  Never use a constant ambient term.  Instead,
  37. *     - For first pixel, cast a whole bunch of rays in different directions
  38. *       from the object intersection point to see what the diffuse illumination
  39. *       really is.  Save this value, after estimating its
  40. *       degree of reusability.  (Method 1)
  41. *     - For second and subsequent pixels,
  42. *         - If there are one or more nearby values already computed,
  43. *           average them and use the result (Method 2), else
  44. *         - Use method 1.
  45. *
  46. *  Implemented by and (c) 1994-6 Jim McElhiney, mcelhiney@acm.org or 71201,1326
  47. *  All standard POV distribution rights granted.  All other rights reserved.
  48. *************************************************************************/
  49.  
  50. #include <string.h>
  51. #include "frame.h"
  52. #include "lighting.h"
  53. #include "vector.h"
  54. #include "povray.h"
  55. #include "optin.h"
  56. #include "povproto.h"
  57. #include "render.h"
  58. #include "texture.h"
  59. #include "octree.h"
  60. #include "radiosit.h"
  61. #include "ray.h"
  62.  
  63.  
  64.  
  65. /*****************************************************************************
  66. * Local preprocessor defines
  67. ******************************************************************************/
  68.  
  69. #define RAD_GRADIENT 1
  70. #define SAW_METHOD 1
  71. /* #define SIGMOID_METHOD 1 */
  72. /* #define SHOW_SAMPLE_SPOTS 1 */   /* try this!  bright spots at sample pts */
  73.  
  74. /*****************************************************************************
  75. * Local typedefs
  76. ******************************************************************************/
  77.  
  78.  
  79.  
  80. /*****************************************************************************
  81. * Local variables
  82. ******************************************************************************/
  83.  
  84.  
  85. long ra_reuse_count = 0;
  86. long ra_gather_count = 0;
  87.  
  88. int Radiosity_Trace_Level = 1;
  89.  
  90. COLOUR Radiosity_Gather_Total;
  91. long Radiosity_Gather_Total_Count;
  92. COLOUR Radiosity_Setting_Total;
  93. long Radiosity_Setting_Total_Count;
  94.  
  95. #ifdef RADSTATS
  96. extern long ot_blockcount;
  97. long ot_seenodecount = 0;
  98. long ot_seeblockcount = 0;
  99. long ot_doblockcount = 0;
  100. long ot_dotokcount = 0;
  101. long ot_lastcount = 0;
  102. long ot_lowerrorcount = 0;
  103. #endif
  104.  
  105.  
  106. OT_NODE *ot_root = NULL;
  107.  
  108. /* This (and all other changing globals) should really be in an execution
  109.  * context structure passed down the execution call tree as a parameter to
  110.  * each function.  This would allow for a multiprocessor/multithreaded version.
  111.  */
  112. FILE *ot_fd = NULL;
  113.  
  114.  
  115. /*****************************************************************************
  116. * Static functions
  117. ******************************************************************************/
  118.  
  119. static long ra_reuse PARAMS((VECTOR IPoint, VECTOR Normal, COLOUR Illuminance));
  120. static long ra_average_near PARAMS((OT_BLOCK *block, void *void_info));
  121. static void ra_gather PARAMS((VECTOR IPoint, VECTOR Normal, COLOUR Illuminance, DBL Weight));
  122. static void VUnpack PARAMS((VECTOR dest_vec, BYTE_XYZ *pack));
  123.  
  124.  
  125.  
  126. /*****************************************************************************
  127. *
  128. * FUNCTION
  129. *
  130. *   Compute_Ambient
  131. *
  132. * INPUT
  133. *   
  134. * OUTPUT
  135. *   
  136. * RETURNS
  137. *   
  138. * AUTHOUR
  139. *
  140. *   Jim McElhiney
  141. *
  142. * DESCRIPTION
  143. *
  144. *   Main entry point for calculated diffuse illumination
  145. *
  146. * CHANGES
  147. *
  148. *   --- 1994 : Creation.
  149. *
  150. ******************************************************************************/
  151.  
  152. int Compute_Ambient(IPoint, S_Normal, Ambient_Colour, Weight)
  153. VECTOR IPoint, S_Normal;
  154. COLOUR Ambient_Colour;   /* the colour to be calculated */
  155. DBL Weight;              /* maximum possible contribution to pixel colour */
  156. {
  157.   int retval, reuse;
  158.   DBL grey, save_bound;
  159.  
  160.   save_bound = opts.Radiosity_Error_Bound;
  161.   if ( Weight < .25 )
  162.   {
  163.     opts.Radiosity_Error_Bound += (.25 - Weight);
  164.   }
  165.   reuse = ra_reuse(IPoint, S_Normal, Ambient_Colour);
  166.   opts.Radiosity_Error_Bound = save_bound;
  167.  
  168.  
  169.   if ( reuse )
  170.   {
  171.     ra_reuse_count++;
  172.     retval = 0;
  173.   }
  174.   else
  175.   {
  176.     ra_gather(IPoint, S_Normal, Ambient_Colour, Weight);
  177.  
  178.     ra_gather_count++;  /* keep a running count */
  179.  
  180.     retval = 1;
  181.   }
  182.  
  183.  
  184.   if ( Radiosity_Trace_Level == 1 )
  185.   {
  186.  
  187.     grey = (Ambient_Colour[RED] + Ambient_Colour[GREEN] + Ambient_Colour[BLUE]) / 3.;
  188.  
  189.     /* note grey spelling:  american options structure with worldbeat calculations! */
  190.     Ambient_Colour[RED]   = opts.Radiosity_Gray * grey + Ambient_Colour[RED]   * (1.-opts.Radiosity_Gray);
  191.     Ambient_Colour[GREEN] = opts.Radiosity_Gray * grey + Ambient_Colour[GREEN] * (1.-opts.Radiosity_Gray);
  192.     Ambient_Colour[BLUE]  = opts.Radiosity_Gray * grey + Ambient_Colour[BLUE]  * (1.-opts.Radiosity_Gray);
  193.  
  194.     /* Scale up by current brightness factor prior to return */
  195.     VScale(Ambient_Colour, Ambient_Colour, opts.Radiosity_Brightness);
  196.   }
  197.  
  198.   return(retval);
  199. }
  200.  
  201.  
  202.  
  203. /*****************************************************************************
  204. *
  205. * FUNCTION
  206. *
  207. *   ra_reuse
  208. *
  209. * INPUT
  210. *
  211. * OUTPUT
  212. *
  213. * RETURNS
  214. *
  215. * AUTHOUR
  216. *
  217. *   Jim McElhiney
  218. *
  219. * DESCRIPTION
  220. *
  221. *   Returns whether or not there were some prestored values close enough to
  222. *   reuse.
  223. *
  224. * CHANGES
  225. *
  226. *   --- 1994 : Creation.
  227. *
  228. ******************************************************************************/
  229.  
  230. static long ra_reuse(IPoint, S_Normal, Illuminance)
  231. VECTOR IPoint, S_Normal;
  232. COLOUR Illuminance;             /* return value */
  233. {
  234.   long i;
  235.   WT_AVG gather;
  236.  
  237.   if (ot_root != NULL)
  238.   {
  239.     Make_Colour(gather.Weights_Times_Illuminances, 0.0, 0.0, 0.0);
  240.  
  241.     gather.Weights = 0.0;
  242.  
  243.     Assign_Vector(gather.P, IPoint);
  244.     Assign_Vector(gather.N, S_Normal);
  245.  
  246.     gather.Weights_Count = 0;
  247.     gather.Good_Count = 0;
  248.     gather.Close_Count = 0;
  249.     gather.Current_Error_Bound = opts.Radiosity_Error_Bound;
  250.  
  251.     for (i = 1; i < Radiosity_Trace_Level; i++)
  252.     {
  253.       gather.Current_Error_Bound *= 1.4;
  254.     }
  255.  
  256.     /*
  257.      * Go through the tree calculating a weighted average of all of the
  258.      * usable points near this one
  259.      */
  260.  
  261.     ot_dist_traverse(ot_root, IPoint, Radiosity_Trace_Level,
  262.                      ra_average_near, (void *)&gather);
  263.  
  264.     /* Did we get any nearby points we could reuse? */
  265.  
  266.     if (gather.Good_Count > 0)
  267.     {
  268.       Make_Colour(gather.Weights_Times_Illuminances, 0.0, 0.0, 0.0);
  269.  
  270.       gather.Weights = 0.0;
  271.  
  272.       for (i = 0; i < gather.Close_Count; i++)
  273.       {
  274.         VAddEq(gather.Weights_Times_Illuminances, gather.Weight_Times_Illuminance[i]);
  275.  
  276.         gather.Weights += gather.Weight[i];
  277.       }
  278.  
  279.       VInverseScale(Illuminance, gather.Weights_Times_Illuminances, gather.Weights);
  280.     }
  281.   }
  282.   else
  283.   {
  284.     gather.Good_Count = 0;      /* No tree, so no reused values */
  285.   }
  286.  
  287.   return(gather.Good_Count);
  288. }
  289.  
  290.  
  291.  
  292. /*****************************************************************************
  293. *
  294. * FUNCTION
  295. *
  296. *   ra_average_near
  297. *
  298. * INPUT
  299. *
  300. * OUTPUT
  301. *
  302. * RETURNS
  303. *
  304. * AUTHOUR
  305. *
  306. *   Jim McElhiney
  307. *   
  308. * DESCRIPTION
  309. *
  310. *   Tree traversal function used by ra_reuse()
  311. *   Calculate the weight of this cached value, taking into account how far
  312. *   it is from our test point, and the difference in surface normal angles.
  313. *
  314. *   Given a node with an old cached value, check to see if it is reusable, and
  315. *   aggregate its info into the weighted average being built during the tree
  316. *   traversal. block contains Point, Normal, Illuminance,
  317. *   Harmonic_Mean_Distance
  318. *
  319. * CHANGES
  320. *
  321. *   --- 1994 : Creation.
  322. *
  323. ******************************************************************************/
  324.  
  325. static long ra_average_near(block, void_info)
  326. OT_BLOCK *block;
  327. void *void_info;
  328. {
  329.   long ind, i;
  330.   WT_AVG *info = (WT_AVG *) void_info;
  331.   VECTOR half, delta, delta_unit;
  332.   COLOUR tc, prediction;
  333.   DBL ri, error_reuse, dir_diff, in_front, dist, weight, square_dist, dr, dg, db;
  334.   DBL error_reuse_rotate, error_reuse_translate, inverse_dist, cos_diff_from_nearest;
  335.   DBL quickcheck_rad;
  336.  
  337.  
  338. #ifdef RADSTATS
  339.   ot_doblockcount++;
  340. #endif
  341.  
  342.   VSub(delta, info->P, block->Point);   /* a = b - c, which is test p minus old pt */
  343.  
  344.   square_dist = VSumSqr(delta);
  345.  
  346.   quickcheck_rad = (DBL)block->Harmonic_Mean_Distance * info->Current_Error_Bound;
  347.  
  348.   /* first we do a tuning test--this func gets called a LOT */
  349.   if (square_dist < quickcheck_rad * quickcheck_rad)
  350.   {
  351.  
  352.     dist = sqrt(square_dist);
  353.     ri = (DBL)block->Harmonic_Mean_Distance;
  354.  
  355.  
  356.     if ( dist > .000001 )
  357.     {
  358.       inverse_dist = 1./dist;
  359.       VScale(delta_unit, delta, inverse_dist);  /* this is a normalization */
  360.  
  361.       /* This block reduces the radius of influence when it points near the nearest
  362.          surface found during sampling. */
  363.       VDot( cos_diff_from_nearest, block->To_Nearest_Surface, delta_unit);
  364.       if ( cos_diff_from_nearest > 0. )
  365.       {
  366.         ri = cos_diff_from_nearest * (DBL)block->Nearest_Distance +
  367.              (1.-cos_diff_from_nearest) * ri;
  368.       }
  369.     }
  370.  
  371.     if (dist < ri * info->Current_Error_Bound)
  372.     {
  373.       VDot(dir_diff, info->N, block->S_Normal);
  374.  
  375.       /* NB error_reuse varies from 0 to 3.82 (1+ 2 root 2) */
  376.  
  377.       error_reuse_translate = dist / ri;
  378.  
  379.       error_reuse_rotate = 2.0 * sqrt(fabs(1.0 - dir_diff));
  380.  
  381.       error_reuse = error_reuse_translate + error_reuse_rotate;
  382.  
  383.       /* is this old point within a reasonable error distance? */
  384.  
  385.       if (error_reuse < info->Current_Error_Bound)
  386.       {
  387. #ifdef RADSTATS
  388.         ot_lowerrorcount++;
  389. #endif
  390.  
  391.         if (dist > 0.000001)
  392.         {
  393.           /*
  394.            * Make sure that the old point is not in front of this point, the
  395.            * old surface might shadow this point and make the result
  396.            * meaningless
  397.            */
  398.           VHalf(half, info->N, block->S_Normal);
  399.           VNormalizeEq(half);          /* needed so check can be to constant */
  400.           VDot(in_front, delta_unit, half);
  401.         }
  402.         else
  403.         {
  404.           in_front = 1.0;
  405.         }
  406.  
  407.         /*
  408.          * Theory:        eliminate the use of old points well in front of our
  409.          * new point we are calculating, but not ones which are just a little
  410.          * tiny bit in front.  This (usually) avoids eliminating points on the
  411.          * same surface by accident.
  412.          */
  413.  
  414.         if (in_front > (-0.05))
  415.         {
  416. #ifdef RADSTATS
  417.           ot_dotokcount++;
  418. #endif
  419.  
  420. #ifdef SIGMOID_METHOD
  421.           weight = error_reuse / info->Current_Error_Bound;  /* 0 < t < 1 */
  422.           weight = (cos(weight * M_PI) + 1.) * .5;           /* 0 < w < 1 */
  423. #endif
  424. #ifdef SAW_METHOD
  425.           weight = 1. - (error_reuse / info->Current_Error_Bound);  /* 0 < t < 1 */
  426. #endif
  427.  
  428.           if ( weight > .001 )
  429.           { /* avoid floating point oddities near zero */
  430.  
  431.             /* This is the block where we use the gradient to improve the prediction */
  432.             dr = delta[X] * block->drdx + delta[Y] * block->drdy + delta[Z] * block->drdz;
  433.             dg = delta[X] * block->dgdx + delta[Y] * block->dgdy + delta[Z] * block->dgdz;
  434.             db = delta[X] * block->dbdx + delta[Y] * block->dbdy + delta[Z] * block->dbdz;
  435. #ifndef RAD_GRADIENT
  436.             dr = dg = db = 0.;
  437. #endif
  438. #if 0
  439.             /* Ensure that the total change in colour is a reasonable magnitude */
  440.             if ( dr > .1 ) dr = .1;  else if ( dr < -.1 ) dr = -.1;
  441.             if ( dg > .1 ) dg = .1;  else if ( dg < -.1 ) dg = -.1;
  442.             if ( db > .1 ) db = .1;  else if ( db < -.1 ) db = -.1;
  443. #endif
  444.             prediction[RED]   = block->Illuminance[RED]   + dr;
  445.             prediction[RED]   = min(max(prediction[RED],  0.), 1.);
  446.             prediction[GREEN] = block->Illuminance[GREEN] + dg;
  447.             prediction[GREEN] = min(max(prediction[GREEN],0.), 1.);
  448.             prediction[BLUE]  = block->Illuminance[BLUE]  + db;
  449.             prediction[BLUE]  = min(max(prediction[BLUE], 0.), 1.);
  450.  
  451. #ifdef SHOW_SAMPLE_SPOTS
  452.             if ( dist < opts.Radiosity_Dist_Max * .015 ) {
  453.               prediction[RED] = prediction[GREEN] = prediction[BLUE] = 3.;
  454.             }
  455. #endif
  456.  
  457.             /* The predicted colour is an extrapolation based on the old value */
  458.             VScale(tc, prediction, weight);
  459.  
  460.             VAddEq(info->Weights_Times_Illuminances, tc);
  461.  
  462.             info->Weights += weight;
  463.  
  464.             info->Weights_Count++;
  465.  
  466.             info->Good_Count++;
  467.  
  468.             ind = -1;
  469.  
  470.             if (info->Close_Count < opts.Radiosity_Nearest_Count)
  471.             {
  472.               ind = info->Close_Count++;
  473.             }
  474.             else
  475.             {
  476.               for (i = 0; i < info->Close_Count; i++)
  477.               {
  478.                 if (dist < info->Distance[i])
  479.                 {
  480.                   ind = i;
  481.  
  482.                   i = opts.Radiosity_Nearest_Count;
  483.                 }
  484.               }
  485.             }
  486.  
  487.             if (ind != -1)
  488.             {
  489.               info->Distance[ind] = dist;
  490.  
  491.               info->Weight[ind] = weight;
  492.  
  493.               VScale(info->Weight_Times_Illuminance[ind], prediction, weight);
  494.  
  495.             }
  496.           }
  497.         }
  498.       }
  499.     }
  500.   }
  501.  
  502.   return(1);
  503. }
  504.  
  505.  
  506.  
  507. /*****************************************************************************
  508. *
  509. * FUNCTION
  510. *
  511. *   ra_gather
  512. *
  513. * INPUT
  514. *   IPoint - a point at which the illumination is needed
  515. *   S_Normal - the surface normal (not perturbed by the current layer) at that point
  516. *   Illuminance - a place to put the return result
  517. *   Weight - the weight of this point in final output, to drive ADC_Bailout
  518. *   
  519. * OUTPUT
  520. *   The average colour of light of objects visible from the specified point.
  521. *   The colour is returned in the Illuminance parameter.
  522. *
  523. *   
  524. * RETURNS
  525. *   
  526. * AUTHOUR
  527. *
  528. *   Jim McElhiney
  529. *   
  530. * DESCRIPTION
  531. *    Gather up the incident light and average it.
  532. *    Return the results in Illuminance, and also cache them for later.
  533. *    Note that last parameter is similar to weight parameter used
  534. *    to control ADC_Bailout as a parameter to Trace(), but it also
  535. *    takes into account that this subsystem calculates only ambient
  536. *    values.  Therefore, coming in at the top level, the value might
  537. *    be 0.3 if the first object hit had an ambient of 0.3, whereas
  538. *    Trace() would have been passed a parameter of 1.0 (since it
  539. *    calculates the whole pixel value).
  540. *
  541. * CHANGES
  542. *
  543. *   --- 1994 : Creation.
  544. *
  545. ******************************************************************************/
  546.  
  547. static void ra_gather(IPoint, S_Normal, Illuminance, Weight)
  548. VECTOR IPoint, S_Normal;
  549. COLOUR Illuminance;
  550. DBL Weight;
  551. {
  552.   extern FRAME Frame;
  553.   long i, hit, Current_Radiosity_Count;
  554.   unsigned long Save_Quality_Flags, Save_Options;
  555.   VECTOR random_vec, direction, n2, n3, up, min_dist_vec;
  556.   int save_Max_Trace_Level;
  557.   DBL Inverse_Distance_Sum, depth, mean_dist, weight, save_min_reuse,
  558.       drdxs, dgdxs, dbdxs, drdys, dgdys, dbdys, drdzs, dgdzs, dbdzs,
  559.       depth_weight_for_this_gradient, dxsquared, dysquared, dzsquared,
  560.       constant_term, deemed_depth, min_dist, reuse_dist_min, to_eye,
  561.       sum_of_inverse_dist, sum_of_dist, average_dist, gradient_count;
  562.   COLOUR Colour_Sums, Temp_Colour;
  563.   RAY New_Ray;
  564.   OT_BLOCK *block;
  565.   OT_ID id;
  566.  
  567.   /*
  568.    * A bit of theory: The goal is to create a set of "random" direction rays
  569.    * so that the probability of close-to-normal versus close-to-tangent rolls
  570.    * off in a cos-theta curve, where theta is the deviation from normal.
  571.    * That is, lots of rays close to normal, and very few close to tangent.
  572.    * You also want to have all of the rays be evenly spread, no matter how
  573.    * many you want to use.  The lookup array has an array of points carefully
  574.    * chosen to meet all of these criteria.
  575.   */
  576.  
  577.  
  578.   /* The number of rays to trace varies with our recursion depth */
  579.  
  580.   Current_Radiosity_Count = opts.Radiosity_Count;
  581.   save_min_reuse = opts.Radiosity_Min_Reuse;
  582.   for ( i=1; i<Radiosity_Trace_Level; i++ )
  583.   {
  584.     Current_Radiosity_Count /= 3;
  585.     opts.Radiosity_Min_Reuse *= 2.;
  586.   }
  587.  
  588.   /* Save some global stuff which we have to change for now */
  589.  
  590.   save_Max_Trace_Level = Max_Trace_Level;
  591.  
  592.  
  593.   /* Since we'll be calculating averages, zero the accumulators */
  594.  
  595.   Make_Colour(Colour_Sums, 0., 0., 0.);
  596.   Inverse_Distance_Sum = 0.;
  597.  
  598.  
  599.   min_dist = BOUND_HUGE;
  600.  
  601.   if ( fabs(fabs(S_Normal[Z])- 1.) < .1 ) {
  602.     /* too close to vertical for comfort, so use cross product with horizon */
  603.     up[X] = 0.; up[Y] = 1.; up[Z] = 0.;
  604.   }
  605.   else
  606.   {
  607.     up[X] = 0.; up[Y] = 0.; up[Z] = 1.;
  608.   }
  609.  
  610.   VCross(n2, S_Normal, up);  VNormalizeEq(n2);
  611.   VCross(n3, S_Normal, n2);  VNormalizeEq(n3);
  612.  
  613.  
  614.   /* Note that this max() forces at least one ray to be shot.
  615.     Otherwise, the loop does nothing, since every call to 
  616.     Trace() just bails out immediately! */
  617.   weight = max(ADC_Bailout, Weight/(DBL)Current_Radiosity_Count);
  618.  
  619.   /* Initialized the accumulators for the integrals which will be come the rad gradient */
  620.   drdxs = dgdxs = dbdxs = drdys = dgdys = dbdys = drdzs = dgdzs = dbdzs = 0.;
  621.   sum_of_inverse_dist = sum_of_dist = gradient_count = 0.;
  622.  
  623.  
  624.   for (i = hit = 0; i < Current_Radiosity_Count; i++)
  625.   {
  626.     /* pick a random direction with the right statistical skew */
  627.  
  628.     VUnpack(random_vec, &rad_samples[i]);
  629.  
  630.     if ( fabs(S_Normal[Z] - 1.) < .001 )         /* pretty well straight Z, folks */
  631.     {
  632.       /* we are within 1/20 degree of pointing in the Z axis. */
  633.       /* use all vectors as is--they're precomputed this way */
  634.       Assign_Vector(direction, random_vec);
  635.     }
  636.     else
  637.     {
  638.       direction[X] = n2[X]*random_vec[X] + n3[X]*random_vec[Y] + S_Normal[X]*random_vec[Z];
  639.       direction[Y] = n2[Y]*random_vec[X] + n3[Y]*random_vec[Y] + S_Normal[Y]*random_vec[Z];
  640.       direction[Z] = n2[Z]*random_vec[X] + n3[Z]*random_vec[Y] + S_Normal[Z]*random_vec[Z];
  641.     }
  642.  
  643.  
  644.     /* Build a ray pointing in the chosen direction */
  645.     Make_Colour(Temp_Colour, 0.0, 0.0, 0.0);
  646.     Initialize_Ray_Containers(&New_Ray);
  647.     Assign_Vector(New_Ray.Initial, IPoint);
  648.     Assign_Vector(New_Ray.Direction, direction);
  649.  
  650.     /* save some flags that must be set to a different value during the trace() */
  651.     Save_Quality_Flags = opts.Quality_Flags;
  652.     Save_Options = opts.Options;
  653.     opts.Radiosity_Quality = 6;
  654.  
  655. #ifdef SAFE_BUT_SLOW
  656.     opts.Quality_Flags = Quality_Values[opts.Radiosity_Quality];
  657. #else
  658.     /* Set up a custom quality level, like Q=5, without area lights, with radiosity */
  659.     opts.Quality_Flags = Q_SHADOW;
  660.     opts.Options &= ~USE_LIGHT_BUFFER;
  661. #endif
  662.  
  663.     /* Go down in recursion, trace the result, and come back up */
  664.  
  665.     Trace_Level++;
  666.     Radiosity_Trace_Level++;
  667.     depth = Trace(&New_Ray, Temp_Colour, weight);
  668.     Radiosity_Trace_Level--;
  669.     Trace_Level--;
  670.  
  671.  
  672.     /* Add into illumination gradient integrals */
  673.  
  674.     deemed_depth = depth;
  675.     if (deemed_depth < opts.Radiosity_Dist_Max * 10.)
  676.     {
  677.       depth_weight_for_this_gradient = 1. / deemed_depth;
  678.       sum_of_inverse_dist += 1. / deemed_depth;
  679.       sum_of_dist += deemed_depth;
  680.       gradient_count++;
  681.  
  682.       dxsquared = direction[X] * direction[X]; if (direction[X] < 0.) dxsquared = -dxsquared;
  683.       dysquared = direction[Y] * direction[Y]; if (direction[Y] < 0.) dysquared = -dysquared;
  684.       dzsquared = direction[Z] * direction[Z]; if (direction[Z] < 0.) dzsquared = -dzsquared;
  685.       drdxs += dxsquared * Temp_Colour[RED]   * depth_weight_for_this_gradient;
  686.       dgdxs += dxsquared * Temp_Colour[GREEN] * depth_weight_for_this_gradient;
  687.       dbdxs += dxsquared * Temp_Colour[BLUE]  * depth_weight_for_this_gradient;
  688.       drdys += dysquared * Temp_Colour[RED]   * depth_weight_for_this_gradient;
  689.       dgdys += dysquared * Temp_Colour[GREEN] * depth_weight_for_this_gradient;
  690.       dbdys += dysquared * Temp_Colour[BLUE]  * depth_weight_for_this_gradient;
  691.       drdzs += dzsquared * Temp_Colour[RED]   * depth_weight_for_this_gradient;
  692.       dgdzs += dzsquared * Temp_Colour[GREEN] * depth_weight_for_this_gradient;
  693.       dbdzs += dzsquared * Temp_Colour[BLUE]  * depth_weight_for_this_gradient;
  694.     }
  695.  
  696.  
  697.     if (depth > opts.Radiosity_Dist_Max)
  698.     {
  699.       depth = opts.Radiosity_Dist_Max;
  700.     }
  701.     else
  702.     {
  703. #ifdef RADSTATS
  704.       hit++;
  705. #endif
  706.     }
  707.  
  708.     if (depth < min_dist)
  709.     {
  710.       min_dist = depth;
  711.       Assign_Vector(min_dist_vec, direction);
  712.     }
  713.  
  714.     opts.Quality_Flags = Save_Quality_Flags;
  715.     opts.Options = Save_Options;
  716.  
  717.     /* Add into total illumination integral */
  718.  
  719.     VAddEq(Colour_Sums, Temp_Colour);
  720.  
  721.     Inverse_Distance_Sum += 1.0 / depth;
  722.   } /* end ray sampling loop */
  723.  
  724.  
  725.   /*
  726.    * Use the accumulated values to calculate the averages needed. The sphere
  727.    * of influence of this primary-method sample point is based on the
  728.    * harmonic mean distance to the points encountered. (An harmonic mean is
  729.    * the inverse of the mean of the inverses).
  730.    */
  731.  
  732.   mean_dist = 1.0 / (Inverse_Distance_Sum / (DBL) Current_Radiosity_Count);
  733.  
  734.   VInverseScale(Illuminance, Colour_Sums, (DBL) Current_Radiosity_Count);
  735.  
  736.   /* Keep a running total of the final Illuminances we calculated */
  737.   if ( Radiosity_Trace_Level == 1.) {
  738.     VAddEq(Radiosity_Gather_Total, Illuminance);
  739.     Radiosity_Gather_Total_Count++;
  740.   }
  741.  
  742.  
  743.   /* We want to cached this block for later reuse.  But,
  744.    * if ground units not big enough, meaning that the value has very
  745.    * limited reuse potential, forget it.
  746.    */
  747.  
  748.   if (mean_dist > (opts.Radiosity_Dist_Max * 0.0001))
  749.   {
  750.  
  751.     /*
  752.      * Theory:  we don't want to calculate a primary method ray loop at every
  753.      * point along the inside edges, so a minimum effectivity is practical.
  754.      * It is expressed as a fraction of the distance to the eyepoint.  1/2%
  755.      * is a good number.  This enhancement was Greg Ward's idea, but the use
  756.      * of % units is my idea.  [JDM]
  757.      */
  758.  
  759.     VDist(to_eye, Frame.Camera->Location, IPoint);
  760.     reuse_dist_min = to_eye * opts.Radiosity_Min_Reuse;
  761.     if (mean_dist < reuse_dist_min)
  762.     {
  763.       mean_dist = reuse_dist_min;
  764.     }
  765.  
  766.  
  767.     /* figure out the block id */
  768.     ot_index_sphere(IPoint, mean_dist * opts.Radiosity_Error_Bound, &id);
  769.  
  770.  
  771. #ifdef RADSTATS
  772.     ot_blockcount++;
  773. #endif
  774.  
  775.     /* After end of ray loop, we've decided that this point is worth storing */
  776.     /* Allocate a block, and fill it with values for reuse in cacheing later */
  777.     block = (OT_BLOCK *)POV_MALLOC(sizeof(OT_BLOCK), "octree block");
  778.     memset(block, 0, sizeof(OT_BLOCK));
  779.  
  780.     /* beta */
  781.     if ( gradient_count > 10.)
  782.     {
  783.       average_dist = sum_of_dist / gradient_count;
  784.       constant_term = 1.00 / (sum_of_inverse_dist * average_dist );
  785.  
  786.       block->drdx = (float)(drdxs * constant_term);
  787.       block->dgdx = (float)(dgdxs * constant_term);
  788.       block->dbdx = (float)(dbdxs * constant_term);
  789.       block->drdy = (float)(drdys * constant_term);
  790.       block->dgdy = (float)(dgdys * constant_term);
  791.       block->dbdy = (float)(dbdys * constant_term);
  792.       block->drdz = (float)(drdzs * constant_term);
  793.       block->dgdz = (float)(dgdzs * constant_term);
  794.       block->dbdz = (float)(dbdzs * constant_term);
  795.     }
  796.  
  797.  
  798.     /* Fill up the values in the octree (ot_) cache block */
  799.  
  800.     Assign_Vector(block->Illuminance, Illuminance);
  801.     Assign_Vector(block->To_Nearest_Surface, min_dist_vec);
  802.     block->Harmonic_Mean_Distance = (float)mean_dist;
  803.     block->Nearest_Distance = (float)min_dist;
  804.     block->Bounce_Depth = (short)Radiosity_Trace_Level;
  805.     Assign_Vector(block->Point, IPoint);
  806.     Assign_Vector(block->S_Normal, S_Normal);
  807.     block->next = NULL;
  808.  
  809.     /* store the info block in the oct tree */
  810.     ot_ins(&ot_root, block, &id);
  811.  
  812.     /* In case the rendering is suspended, save the cache tree values to a file */
  813.     if ( opts.Radiosity_File_SaveWhileRendering && (ot_fd != NULL) ) {
  814.       ot_write_block(block, ot_fd);
  815.     }
  816.   }
  817.  
  818.   /* Put things back where they were in recursion depth */
  819.   Max_Trace_Level = save_Max_Trace_Level;
  820.   opts.Radiosity_Min_Reuse = save_min_reuse;
  821.  
  822. }
  823.  
  824.  
  825. /*****************************************************************************
  826. *
  827. * FUNCTION  VUnpack()  -  Unpacks "pack_vec" into "dest_vec" and normalizes it.
  828. *
  829. * INPUT
  830. *
  831. * OUTPUT
  832. *
  833. * RETURNS   Nothing
  834. *
  835. * AUTHOUR   Jim McElhiney
  836. *
  837. * DESCRIPTION
  838. *
  839. *  The precomputed radiosity rays are packed into a lookup array with one byte
  840. *  for each of dx, dy, and dz.  dx and dy are scaled from the range (-1. to 1.),
  841. *  and dz is scaled from the range (0. to 1.), and both are stored in the range
  842. *  0 to 255.
  843. *
  844. *  The reason for this function is that it saves a bit of memory.  There are 2000
  845. *  entries in the table, and packing them saves 21 bytes each, or 42KB.
  846. *
  847. * CHANGES
  848. *
  849. *   --- Jan 1996 : Creation.
  850. *
  851. ******************************************************************************/
  852. static void VUnpack(dest_vec, pack_vec)
  853. VECTOR dest_vec;
  854. BYTE_XYZ * pack_vec;
  855. {
  856.   dest_vec[X] = ((double)pack_vec->x * (1./ 255.))*2.-1.;
  857.   dest_vec[Y] = ((double)pack_vec->y * (1./ 255.))*2.-1.;
  858.   dest_vec[Z] = ((double)pack_vec->z * (1./ 255.));
  859.  
  860.   VNormalizeEq(dest_vec);   /* already good to about 1%, but we can do better */
  861. }
  862.  
  863.  
  864. /*****************************************************************************
  865. *
  866. * FUNCTION  Initialize_Radiosity_Code
  867. *
  868. * INPUT     Nothing.
  869. *
  870. * OUTPUT    Sets various global states used by radiosity.  Notably,
  871. *           ot_fd - the file identifier of the file used to save radiosity values
  872. *
  873. * RETURNS   1 for Success, 0 for failure  (e.g., could not open cache file)
  874. *
  875. * AUTHOUR   Jim McElhiney
  876. *
  877. * DESCRIPTION
  878. *
  879. * CHANGES
  880. *
  881. *   --- Jan 1996 : Creation.
  882. *
  883. ******************************************************************************/
  884. long
  885. Initialize_Radiosity_Code()
  886. {
  887.   long retval, used_existing_file;
  888.   FILE *fd;
  889.   char *modes, rad_cache_filename[256];
  890.  
  891.   retval = 1;                   /* assume the best */
  892.  
  893.   if ( opts.Options & RADIOSITY )
  894.   {
  895.     opts.Radiosity_Preview_Done = 0;
  896.  
  897.     ra_gather_count  = 0;
  898.     ra_reuse_count   = 0;
  899.  
  900.  
  901.     if ( opts.Radiosity_Dist_Max == 0. )
  902.     {
  903.       /* User hasn't picked a radiosity dist max, so pick one automatically. */
  904.       VDist(opts.Radiosity_Dist_Max, Frame.Camera->Location,
  905.                                      Frame.Camera->Look_At);
  906.       opts.Radiosity_Dist_Max *= 0.2;
  907.     }
  908.  
  909. #ifdef RADSTATS
  910.     ot_seenodecount  = 0;
  911.     ot_seeblockcount = 0;
  912.     ot_doblockcount  = 0;
  913.     ot_dotokcount    = 0;
  914.     ot_lowerrorcount = 0;
  915.     ot_lastcount     = 0;
  916. #endif
  917.  
  918.     if ( ot_fd != NULL )    /* if already open for some unknown reason, close it */
  919.     {
  920.       fclose(ot_fd);
  921.       ot_fd = 0;
  922.     }
  923.  
  924.     /* build the file name for the radiosity cache file */
  925.     strcpy(rad_cache_filename, opts.Scene_Name);
  926.     strcat(rad_cache_filename, RADIOSITY_CACHE_EXTENSION);
  927.  
  928.     used_existing_file = 0;
  929.     if ( ((opts.Options & CONTINUE_TRACE) && opts.Radiosity_File_ReadOnContinue)  ||
  930.          opts.Radiosity_File_AlwaysReadAtStart )
  931.     {
  932.       fd = fopen(rad_cache_filename, READ_FILE_STRING);   /* "myname.rca" */
  933.       if ( fd != NULL) {
  934.         used_existing_file = ot_read_file(fd);
  935.         retval &= used_existing_file;
  936.         fclose(fd);
  937.       }
  938.     }
  939.     else
  940.     {
  941.       DELETE_FILE(rad_cache_filename);  /* default case, force a clean start */
  942.     }
  943.  
  944.  
  945.     if ( opts.Radiosity_File_SaveWhileRendering )
  946.     {
  947.       /* If we are writing a file, but not using what's there, we truncate,
  948.          since we conclude that what is there is bad.
  949.          But, if we are also using what's there, then it must be good, so
  950.          we just append to it.
  951.       */
  952.       modes = used_existing_file ? APPEND_FILE_STRING : WRITE_FILE_STRING;
  953.       ot_fd = fopen(rad_cache_filename, modes);
  954.       retval &= (ot_fd != NULL);
  955.     }
  956.   }
  957.  
  958.   return retval;
  959. }
  960.  
  961.  
  962. /*****************************************************************************
  963. *
  964. * FUNCTION  Deinitialize_Radiosity_Code()
  965. *
  966. * INPUT     Nothing.
  967. *
  968. * OUTPUT    Sets various global states used by radiosity.  Notably,
  969. *           ot_fd - the file identifier of the file used to save radiosity values
  970. *
  971. * RETURNS   1 for total success, 0 otherwise (e.g., could not save cache tree)
  972. *
  973. * AUTHOUR   Jim McElhiney
  974. *
  975. * DESCRIPTION
  976. *   Wrap up and free any radiosity-specific features.
  977. *   Note that this function is safe to call even if radiosity was not on.
  978. *
  979. * CHANGES
  980. *
  981. *   --- Jan 1996 : Creation.
  982. *
  983. ******************************************************************************/
  984. long
  985. Deinitialize_Radiosity_Code()
  986. {
  987.   long retval;
  988.   char rad_cache_filename[256];
  989.   FILE *fd;
  990.  
  991.   retval = 1;                    /* assume the best */
  992.  
  993.   /* if the global file identifier is set, close it */
  994.   if ( ot_fd != NULL ) {
  995.     fclose(ot_fd);
  996.     ot_fd = NULL;
  997.   }
  998.  
  999.  
  1000.   /* build the file name for the radiosity cache file */
  1001.   strcpy(rad_cache_filename, opts.Scene_Name);
  1002.   strcat(rad_cache_filename, RADIOSITY_CACHE_EXTENSION);
  1003.  
  1004.  
  1005.   /* If user has not asked us to save the radiosity cache file, delete it */
  1006.   if ( opts.Radiosity_File_SaveWhileRendering  &&
  1007.       !(opts.Radiosity_File_KeepAlways || (Stop_Flag && opts.Radiosity_File_KeepOnAbort) ) )
  1008.   {
  1009.     DELETE_FILE(rad_cache_filename);
  1010.   }
  1011.  
  1012.   /* after-the-fact version.  This is an alternative to putting a call to 
  1013.      ot_write_node after the call to ot_ins in ra_gather().
  1014.      The on-the-fly version (all of the code which uses ot_fd) is superior
  1015.      in that you will get partial results if you restart your rendering
  1016.      with a different resolution or camera angle.  This version is superior
  1017.      in that your rendering goes a lot quicker.
  1018.   */
  1019.  
  1020.   if (!(opts.Radiosity_File_KeepAlways || (Stop_Flag && opts.Radiosity_File_KeepOnAbort)) &&  
  1021.       !opts.Radiosity_File_SaveWhileRendering )
  1022.   {
  1023.     fd = fopen(rad_cache_filename, WRITE_FILE_STRING);
  1024.  
  1025.     if ( fd != NULL ) {
  1026.       retval &= ot_save_tree(ot_root, fd);
  1027.  
  1028.       fclose(fd);
  1029.     }
  1030.     else
  1031.     {
  1032.       retval = 0;
  1033.     }
  1034.   }
  1035.  
  1036.  
  1037.   /* Note that multiframe animations should call this free function if they have
  1038.      moving objects and want correct results.
  1039.      They should NOT call this function if they have no moving objects (like
  1040.      fly-throughs) and want speed
  1041.   */
  1042.   if ( ot_root != NULL ) {
  1043.     retval &= ot_free_tree(&ot_root);   /* this zeroes the root pointer */
  1044.   }
  1045.  
  1046.   return retval;
  1047. }
  1048.